Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 
http://www.biomedcentral.eom/1 471 -21 05/1 3/S5/S8 



RESEARCH Open Access 



Exploring behaviors of stochastic differential 
equation models of biological systems using 
change of measures 

Sumit Kumar Jha 1 ", Christopher James Langmead 2,3 " 

From First IEEE International Conference on Computational Advances in Bio and medical Sciences (ICCABS 
2011) 

Orlando, FL, USA. 3-5 February 201 1 



Abstract 

Stochastic Differential Equations (SDE) are often used to model the stochastic dynamics of biological systems. 
Unfortunately, rare but biologically interesting behaviors (e.g., oncogenesis) can be difficult to observe in stochastic 
models. Consequently, the analysis of behaviors of SDE models using numerical simulations can be challenging. 
We introduce a method for solving the following problem: given a SDE model and a high-level behavioral 
specification about the dynamics of the model, algorithmically decide whether the model satisfies the specification. 
While there are a number of techniques for addressing this problem for discrete-state stochastic models, the 
analysis of SDE and other continuous-state models has received less attention. Our proposed solution uses a 
combination of Bayesian sequential hypothesis testing, non-identically distributed samples, and Girsanov's theorem 
for change of measures to examine rare behaviors. We use our algorithm to analyze two SDE models of tumor 
dynamics. Our use of non-identically distributed samples sampling contributes to the state of the art in statistical 
verification and model checking of stochastic models by providing an effective means for exposing rare events in 
SDEs, while retaining the ability to compute bounds on the probability that those events occur. 



Background 

The dynamics of biological systems are largely driven by 
stochastic processes and subject to random external per- 
turbations. The consequences of such random processes 
are often investigated through the development and ana- 
lysis of stochastic models (e.g., [1-4]). Unfortunately, the 
validation and analysis of stochastic models can be very 
challenging [5,6], especially when the model is intended 
to investigate rare, but biologically significant behaviors 
(e.g., oncogenesis). The goal of this paper is to introduce 
an algorithm for examining such rare behaviors in Sto- 
chastic Differential Equation (SDE) models. The algo- 
rithm takes as input the SDE model, M. > and a high- 
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level description of a dynamical behavior, cp (e.g., a for- 
mula in temporal logic). It then computes a statistically 
rigorous bound on the probability that the given model 
exhibits the stated behavior using a combination of 
biased sampling and Bayesian Statistical Model Check- 
ing [7,8]. 

Existing methods for validating and analyzing stochas- 
tic models often require extensive Monte Carlo sam- 
pling of independent trajectories to verify that the 
model is consistent with known data, and to character- 
ize the model's expected behavior under various initial 
conditions. Sampling strategies are either unbiased or 
biased. Unbiased sampling strategies draw trajectories 
according to the probability distribution implied by the 
model, and are thus not well-suited to investigating rare 
behaviors. For example, if the actual probability that the 
model will exhibit a given behavior is 10~ 10 , then the 
expected number of samples need to see such behaviors 
is about 10 10 (See Figure 1). Biased sampling strategies, 
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Figure 1 Observing rare behaviors in i.i.d. sampling is challenging. A toy model with a one low-probability state. An unbiased sampling 
algorithm may require billions of samples in order to observe the 'bad' state. Statistical algorithms based on i.i.d. algorithms are not suitable for 
analyzing such models with rare interesting behaviors. 



in contrast, can be used to increase the probability of 
observing rare events, at the expense of distorting the 
underlying probability distribution/measure. If the 
change of measure is well-characterized (as in impor- 
tance sampling), these distortions can be corrected 
when computing properties of the distribution. 

Our method uses a combination of biased sampling 
and sequential hypothesis testing [9,10] to estimate the 
probability that the model satisfies the property. Briefly, 
the algorithm randomly perturbs the computational 
model prior to generating each sample in order to expose 
rare but interesting behaviors. These perturbations cause 
a change of measure. We note that in the context of 
SDEs, a change of measure is itself a stochastic process 
and so importance sampling, which assumes that the 
change of measure with respect to the biased distribution 
is known, cannot be used. Our technique does not 
require us to know the exact magnitudes of the changes 
of measure, nor the Radon-Nikodym derivatives [11]. 
Instead, it ensures that that the geometric average of 
these derivatives is bounded. This is a much weaker 
assumption than is required by importance sampling, but 
we will show that it is sufficient for the purposes of 
obtaining bounds via sequential hypothesis testing. 

Related work 

Our method performs statistical model checking using 
hypothesis testing [12,13], which has been used pre- 
viously to analyze in a variety of domains (e.g., 
[14-19,19]), including computational biology (e.g., 



[7,8,20-22]). There are two key differences between the 
method in this paper and existing methods. First, exist- 
ing methods generate independently and identically dis- 
tributed (i.i.d.) samples. Our method, in contrast, 
generates independent but non-identically distributed 
(non- i.i.d.) samples. It does so to expose rare behaviors. 
Second, whereas most existing methods for statistical 
model checking use classical statistics, our method 
employs Bayesian statistics [23,24], We have previously 
shown that Bayesian statistics confers advantages in the 
context of statistical verification in terms of efficiency 
and flexibility [7,8,19,25]. 

Methods 

Our method draws on concepts from several different 
fields. We begin by briefly surveying the semantics of 
stochastic differential equations, a language for formally 
specifying dynamic behaviors, Girsanov's theorem on 
change of measures, and results on consistency and con- 
centration of Bayesian posteriors. 

Stochastic differential equation models 

A stochastic differential equation (SDE) [26,27] is a dif- 
ferential equation in which some of the terms evolve 
according to Brownian Motion [28]. A typical SDE is of 
the following form: 

dX = b{t, X t )dt + v[t, X t )dW t 

where X is a system variable, b is a Riemann integr- 
able function, v is an Ito integrable function, and W is 
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Brownian Motion. The Brownian Motion W is a contin- 
uous-time stochastic process satisfying the following 
three conditions: 

1. Wo = 0 

2. W t is continuous {almost surely), 

3. W t has independent normally distributed 
increments: 

• W t 'W s and W f - W s > are independent if 0 < s 

<t <s <t\ 

♦ W t - W s ~Af{0, t-s), where J\f{0, t - s) 
denotes the normal distribution with mean 0 and 
variance t - s. Note that the symbol ~ is used to 
indicate "is distributed as". 

Consider the time between 0 and t as divided into m 
discrete steps 0, t lt t 2 ... t m = t. The solution to a sto- 
chastic differential equation is the limit of the following 
discrete difference equation, as m goes to infinity: 

Xt k+ i ~X tk = b(tk, X tk ) (tfe+i - tk) 

+ v{t k ,X tk ){W tk+1 -W tk ) 

In what follows, M. will refer to a system of stochas- 
tic differential equations. We note that a system of sto- 
chastic differential equations comes equipped with an 
inherent probability space and a natural probability 
measure Our algorithm repeatedly and randomly per- 
turbs the probability measure of the Brownian motion 
in the model J\A which, in turn, changes the underlying 
measure in an effort to expose rare behaviors. These 
changes can be characterized using Girsanovs Theorem. 
Girsanov's theorem for perturbing stochastic differential 
equation models 

Given a process {6 t | 0 < t < T) satisfying the Novikov 
condition [29], such as an SDE, the following exponen- 
tial martingale Z t defines the change from measure P to 
new measure p : 

Z t = exp{- [ O u dW u - [ 0^du/2) 
Jo Jo 

Here, Z t is the Radon-Nikodym derivative of p with 
respect to P for t <T. The Brownian motion w t under 
p is given by: W t = W t + f^O u du. The non-stochastic 
component of the stochastic differential equation is not 
affected by change of measures. Thus, a change of mea- 
sure for SDEs is a stochastic process (unlike importance 
sampling for explicit probability distributions). 

Specifying dynamic behaviors 

Next, we define a formalism for encoding high-level 
behavioral specification that our algorithm will test 
against M ♦ 



Definition 1 (Adapted Finitely Monitorable). Let a be 
a finite-length trace from the stochastic differential 
equation M . A specification cp is said to be adapted 
finitely monitorable (AFM) if it is possible to decide 
whether a satisfies (p, denoted (p \= (p. 

Certain AFM specifications can be expressed as for- 
mulas in Bounded Linear Temporal Logic (BLTL) 
[30-32]. Informally, BLTL formulas can capture the 
ordering of events. 

Definition 2 (Probabilistic Adapted Finitely Monitor- 
able). A specification (p is said to be probabilistic 
adapted finitely monitorable (PAFM) if it is possible to 
(deterministically or probabilistically) decide whether 
M satisfies cp with probability at least 6, denoted 
M\=P> e {4>). 

Some common examples of PAFM specifications 
include Probabilistic Bounded Linear Temporal Logic 
(PBLTL) (e.g., see [19]) and Continuous Specification 
Logic. Note that temporal logic is only one means for 
constructing specifications; other formalisms can also be 
used, like Statecharts [33,34]. 
Semantics of bounded linear temporal logic (BLTL) 
We define the semantics of BLTL with respect to the 
paths of M ♦ Let a = (s 0 , A 0 ), (s lt A!),... be a sampled 
execution of the model along states s 0 , s lf ... with dura- 
tions A 0 , Ax, .. . e EL We denote the path starting at 
state i by a 1 (in particular, a 0 denotes the original 
execution a). The value of the state variable x in a at 
the state i is denoted by V (or, i, The semantics of 
BLTL is defined as follows: 

1= x ~ v if and only if Via, k, x) ~ v, where v e 
R and ~ e {>, <, =}. 

2. o* i= <pi yq> 2 if and only if a* i= (pi or o* i= (p 2 . 
\= q>i A (p 2 if and only if o* i= q> x and o* 1= (p 2 ; 

4. a* \= -i^>! if and only if a* 1= q>i does not hold. 

5. a* i= <piU*<p2 if an d on ly if there exists i e N such 
that: (a) 0 < Z 0 < / </ A^ +/ < t; (b) a k+i \= (p 2 ; and (c) for 
each 0 < ; < /, d <+j \= q>i, 

Statistical model validation 

Our algorithm performs statistical model checking using 
Bayesian sequential hypothesis testing [23] on non-LLd. 
samples. The statistical model checking problem is to 
decide whether a model M satisfies a probabilistic 
adapted finitely monitorable formula cp with probability 
at least 0. That is, whether M N P>#(0) where 6 e (0, 
1). 

Sequential hypothesis testing 

Let p be the (unknown but fixed) probability of the 
model satisfying cp. We can now re-state the statistical 
model checking problem as deciding between the two 
composite hypotheses: H 0 : p ^ 6 and Hi : p < 6, Here, 
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the null hypothesis H 0 indicates that M satisfies the 
AFM formula cp with probability at least 6, while the 
alternate hypothesis H x denotes that M satisfies the 
AFM formula (p with probability less than 6. 

Definition 3 (Type I and II errors). A Type I error is 
an error where the hypothesis test asserts that the null 
hypothesis H 0 is false, when in fact H 0 is true. Conver- 
sely, a Type II error is an error where the hypothesis 
test asserts that the null hypothesis H 0 is true, when in 
fact H 0 is false. 

The basic idea behind any statistical model checking 
algorithm based on sequential hypothesis testing is to 
iteratively sample traces from the stochastic process. 
Each trace is then evaluated with a trace verifier [32], 
which determines whether the trace satisfies the specifi- 
cation i.e. Gi i= (p. This is always feasible because the spe- 
cifications used are adapted and finitely monitorable. 
Two accumulators keep track of the total number of 
traces sampled, and the number of satisfying traces, 
respectively. The procedure continues until there is 
enough information to reject either H 0 or H\. 
Bayesian sequential hypothesis testing 
Recall that for any finite trace G t of the system and an 
Adapted Finitely Monitorable (AFM) formula (p, we can 
decide whether G t satisfies (p. Therefore, we can define a 
random variable X t denoting the outcome of G t i= (p. 
Thus, Xi is a Bernoulli random variable with probability 
mass function f(xi\p) = p Xi (l — p) x ~ Xi > where x t = 1 if 
and only if G t \= q (p } otherwise x t = 0. 

Bayesian statistics requires that prior probability distri- 
butions be specified for the unknown quantity, which is 
p in our case. Thus we will model the unknown quan- 
tity as a random variable u with prior density g(u). The 
prior probability distribution is usually based on our 
previous experiences and beliefs about the system. Non- 
informative or objective prior probability distribution 
[35] can be used when nothing is known about the 
probability of the system satisfying the AFM formula. 

Suppose we have a sequence of independent random 
variables Xi,..., X n defined as above, and let d = (x lf ..., 
x n ) denote a sample of those variables. 

Definition 4. The Bayes factor B of sample d and 

hypotheses H 0 and H 1 is o = y 

The Bayes factor may be used as a measure of relative 
confidence in H 0 vs. H 1} as proposed by Jeffreys [36]. 
The Bayes factor is the ratio of two likelihoods: 

B = f Q 1 f(xi\u)---f(x n \u)-g(u)du 
J e Q f(x l \u)---f(x n \u) -g(u)du 

We note that the Bayes factor depends on both the 
data d and on the prior g, so it may be considered a 



measure of confidence in H 0 vs. H Y provided by the data 

x w and "weighted" by the prior g. 
Non-U.d. Bayesian sequential hypothesis testing 

Traditional methods for hypothesis testing, including 
those outlined in the previous two subsection*s, assume 
that the samples are drawn Ltd.. In this section" we show 
that non-i.i.d. samples can also be used, provided that cer- 
tain conditions hold. In particular, if one can bound the 
change in measure associated with the non-identical sam- 
pling, one can can also bound the Type I and Type II 
errors under a change of measure. Our algorithm bounds 
the change of measure, and thus the error. 

We begin by reviewing some fundamental concepts 
from Bayesian statistics including KL divergence, KL 
support, affinity, and ^-separation, and then restate an 
important result on the concentration of Bayesian pos- 
teriors [35,37]. 

Definition 5 (Kullback-Leibler (KL) Divergence). Given 
a parameterized family of probability distributions {/#}, 
the Kullback-Leibler (KL) divergence K(6 0 , 6) between 
the distributions corresponding to two parameters 6 and 
6 0 is: K (Go, 0) = Eq 0 \felfe Q ] . Note that Eq q is the expecta- 
tion computed under the probability measure fe 0 . 

Definition 6 (KL Neighborhood). Given a parameter- 
ized family of probability distributions {/#}, the KL 
neighborhood K e (6 0 ) of a parameter value 6 0 is given by 
the set {9 : K (6 0 , 6) < e}. 

Definition 7 (KL Support). A point 9 0 is said to be in 
the KL support of a prior II if and only if for all s >0, II 
(K e (0 0 )) >0. 

Definition 8 (Affinity). The affinity Aff(/i g) between 
any two densities is defined as Aff(f, g) = J y[jgA\i . 

Definition 9 (Strong Separation). Let i c [0, 1] and 
S >0. The set A and the point 9 0 are said to be strongly 
^-separated if and only if for any proper probability dis- 
tribution v on A, Afi(fo 0 , f A fe(xi)v(dO)) < 8 . 

Given these definitions, it can be shown that the Baye- 
sian posterior concentrates exponentially under certain 
technical conditions [35,37]. 
Bounding errors under a change of measure 
Next, we develop the machinery needed to compute 
bounds on the Type-I/Type-II errors for a testing strat- 
egy based on non-LLd. samples. 

A stochastic differential equation model M is natu- 
rally associated with a probability measure Our non-/. 
Id. sampling strategy can be thought of as the assign- 
ment of a set of probability measures ft lt {t 2 ><- to M. . 
Each unique sample G t is associated with an implied 
probability measure fa and is generated from M under 
(Ai in an Lid. manner. Our proofs require that all the 
implied probability measures are equivalent. That is, an 
event is possible (resp. impossible) under a probability 
measure if and only if it is possible (resp. impossible) 
under the original probability measure ft. 
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We use the following result regarding change of mea- 
sures. Suppose a given behavior, say (p, holds on the ori- 
ginal model with an (unknown) probability p. 



P(p<0 0 \X i ) = 



fo° PnjXj\u) g(u) du 
foPv,{Xi\u) g(u) du 



Here, X t is a Bernoulli random variable denoting the 



event that i sample satisfies the given behavior cp. Note 
that the X t s must be independent of one another. Now, 
we can rewrite the above expression as: 



P{p <0 0 \Xi) 



Pn{Xj\u) 
Pm(Xi\u) 



Note that the term p^X^u) denotes the probability of 
observing the event X t under the modified probability 
measure ^ if the unknown probability p were u. In 
order to ensure the independence assumption, the new 
probability measure^ /^^g ^chosen independently of one 
another. The ratio Ml 1 is the implied Radon-Niko- 
dym derivative for ute dhange of measure between two 
equivalent probability measures. Suppose, the testing 
strategy has made n observations X lt X 2 >... X n . Then, 



P(p<0 0 \X) 



Jo n I rv^ fc^M gfa) du 

i=l P^AilU) 



foil 

i=l 



on, fY . A Pm{Xi\u) g{u) du 



A sampling algorithm can compute p^X^u) by 
drawing independent samples from a stochastic differen- 
tial equation model under the new "modified" probabil- 
ity measure. We note that it is not easy to compute the 
change of measure algebraically or numerically. 

However, our algorithm does not need to compute this 
quantity explicitly. It simply establishes bounds on it. 

Consider the following expression that is computable 
without knowing the implied Radon-Nikodym derivative 
or change of measure explicitly. 



Q(p<e 0 \Xi) 



fo° Pmi X i\ u ) %{ u ) du 



foPm( x i\u) g(u) du 
Now, we can rewrite the above expression as: 



Q(p <o 0 \Xi) 



Jo P^ x i\u) 



Pti{Xi\u) g{u) du 



/o 1 



P^{Xj\u) 
P^ x i\u) 



p^{Xi\u) g{u) du 



Our result will exploit the fact that we do not allow 
our testing or sampling procedures to have arbitrary 
implied Radon-Nikodym derivatives. This is reasonable 
as no statistical guarantees should be available for an 
intelligently designed but adversarial test procedure that 
(say) tries to avoid sampling from the given behavior. 
Suppose that the implied Radon-Nikodym derivative 
always lies between a constant c and another constant 
l/c. That is, the change of measure does not distort the 
probabilities of observable events by more than a factor 
of c. Then, we observe that: 

Q{p<0 0 \X l )<I^ MXilu)g{u)dU 

rl I 



Jo -Pv( X i\ U )g( U ) du 

= c 2 p{p<e Q \x i ) 



Furthermore, 



Q(p <o 0 \Xi) > P{ p< e 0 \Xi). 

c 1 

Thus, by allowing the sampling algorithm to change 
measures by at most c, we have changed the posterior 
probability of observing a behavior by at most c 2 . 

Example: Suppose, the testing strategy has made n 
observations X lt X 2 >... X n . Then, 



Q{p<0 0 \X) = 



S:°Y\ P f^P^\u)g{u)du 

i=l PfiV A i\ u ) 



\ Pm{Xi\u) 



i=l \>\i 



Pn(Xi\u) g{u) du 



fo 0 c n UiP^Xi\u))g{u)du 



i=i 



Similarly, 



Q(P <0 0 \X) 



~ i 1 " 

Jo -^UiPn{Xi\u))g{u)du 

= c 2n P('p <0 0 \X U X 2l ...X n ) 



, fo°U{P,m^)s(u)du 

1 1=1 

r 2n i n 

" loU{p,mu))g(u)du 

i=i 



r 2n 



P{p<0 0 \X l , ...X n ) 



Termination conditions for non-U.d. sampling 

Traditional (i.e., Ltd.) Bayesian Sequential Hypothesis 
Testing is guaranteed to terminate. That is, only a finite 
number of samples are required before the test selects 
one of the hypotheses. We now consider the conditions 
under which a Bayesian Sequential Hypothesis Testing 
based procedure using non-i.i.d. samples will terminate. 
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To do this, we first need to show that the posterior 
probability distribution will concentrate on a particular 
value as we see more an more samples from the model 

To consider the conditions under which our algorithm 
will terminate after observing n samples, note that the 
factor introduced due to the change of measure c 2n can 
outweigh the gain made by the concentration of the 
probability measure e' nh . This is not surprising because 
our construction thus far does not force the test not to 
bias against a sample in an intelligent way. That is, a 
maliciously designed testing procedure could simply 
avoid the error prone regions of the design. To address 
this, we define the notion of a fair testing strategy that 
does not engage in such malicious sampling. 

Definition 10. A testing strategy is 77-fair (77 > 1) if 
and only if the geometric average of the implied Radon- 
Nikodym derivatives over a number of samples is within 
a constant factor 77 of unity, i.e., 



n 



Pn{Xi\u) 
Pn(Xi\u) 



Note that a fair test strategy does not need to sample 
from the underlying distribution in an Ltd. manner. 
However, it must guarantee that the probability of 
observing the given behavior in a large number of 
observations is not altered substantially by the non-LLd. 
sampling. Intuitively, we want to make sure that we bias 
for each sample as many times as we bias against it. 
Our main result shows that such a long term neutrality 
is sufficient to generate statistical guarantees on an 
otherwise non-lid. testing procedure. 

Definition 11. An 77-fair test is said to be eventually 
fair if and only if 1 < if <e b , where b is the constant in 
the exponential posterior concentration theorem. 

The notion of a eventually fair test corresponds to a 
testing strategy that is not malicious or adversarial, and 
is making an honest attempt to sample from all the 
events in the long run. 

Algorithm 

Finally, we present our Statistical Verification algo- 
rithm (See Figure 2) in terms of a generic non-LLd. 
testing procedure sampling with random "implied" 
change of measures. Our algorithm is relatively simple 
and generalizes our previous Bayesian Statistical verifi- 
cation algorithm [8] to non-LLd. samples using change 
of measures. The algorithm draws non-LLd. samples 
from the stochastic differential equation under ran- 
domly chosen probability measures. The algorithm 
ensures that the implied change of measure is bounded 
so as to make the testing approach fair. The variable n 
denotes the number of samples obtained so far and x 



denoted the number of samples that satisfy the AFM 
specification (p. Based upon the samples observed, we 
compute the Bayes Factor under the new probability 
measures. We know that the Bayes Factor so computed 
is within a factor of the original Bayes Factor under 
the natural probability measure. Hence, the algorithm 
divides the Bayes Factor by the factor rf n if the Bayes 
Factor is larger than one. If the Bayes Factor is less 
than one, the algorithm multiplies the Bayes Factor by 
the factor if n . 

Results and discussion 

We applied our algorithm to two SDE models of tumor 
dynamics from the literature. The first model is a single 
dimensional stochastic differential equation for the 
influence of chemotherapy on cancer cells, and the sec- 
ond model is a pair of SDEs that describe an immuno- 
genic tumor. 

Lefever and Garay model 

Lefever and Garay [38] studied the growth of tumors 
under immune surveillance and chemotherapy using the 
following stochastic differential equation: 



dx f x fix 2 
— = 1\)X ( 1 ) - 

+ x(l-|)W t 



+ X(l ) AqCOS {cot) 

K 



Here, x is the amount of tumor cells, A 0 cos(cot) 
denotes the influence of a periodic chemotherapy treat- 
ment, r 0 is the linear per capita birth rate of cancer 
cells, K is the carrying capacity of the environment, and 
P represents the influence of the immune cells. Note 
that W t is the standard Brownian Motion, and default 
model parameters were those used in [38].. 

We demonstrate our algorithm on a simple property 
of the model. Namely, starting with a tumor consisting 
of a billion cells, is there at least a 1% chance that the 
tumor could increase to one hundred billion cells under 
under immune surveillance and chemotherapy. The fol- 
lowing BLTL specification captures the behavioral speci- 
fication: 

Pr> 0 .oi(F 10 (x> 10 11 )) 

Figure 3 contrasts the number of samples needed to 
decide whether the model satisfies the property using Li. 
d. and non-LLd.. As expected, the number of samples 
required increases linearly in the logarithm of the Bayes 
factor regardless of whether i.i.d. or non-LLd. sampling 
is used. However, non-LLd. sampling always requires 
fewer samples than i.i.d. sampling. Moreover, the differ- 
ence between the number of samples increases with the 
Bayes factor. That is, the lines are diverging. 
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Require: stochastic differential equation Model M. AFM Specification 0, Bound on implied Radon- Nikodym derivatives r). 
Ensure: Probability of Type I/II error is less than S. 

Initialize ra=l; /*n is number of samples observed*/ 
Initialize .r=0: /*./• is number of successful samples*/ 
B <- 1 

while ( B > \ and B < 6 ) do 

Choose a random implied change of probability measure c, using an ^-eventually fair testing strategy. 
Draw sample a n from the stochastic differential equation under the new probability measure. 

if a n (= 0 then 

X n <- 1 ; x 4- x + 1 
else 

X n <- 0 : 
end if 

ri 



r 



B 



i 



nPin(Xi\u)g{u) du 



if B > 1 then 

B 4- -\-B 

else 

B 4- if"B 
end if 

n 4- n + 1 ; 
end while 

Figure 2 Non-i.i.d. Statistical Verification Algorithm. The figure illustrated the non-i.i.d. Bayesian model validation algorithm. The algorithm 
builds upon Girsanov's theorem on change of measure and Bayesian model validation. 
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We note that there are circumstances when our algo- 
rithm may require more samples than one based on Li. 
d. sampling. This will happen when the property being 
tested has relatively high probability. For example, we 
tested the property that the probability of eradicating 
the tumor is at most 1%. The i.i.d. algorithm required 
1374 samples to deny this possibility while the non-i.i.d. 
algorithm required 1526 samples. Thus, our algorithm is 
best used to examine rare behaviors. 

Nonlinear immunogenic tumor model 

The second model we analyze studies immunogenic 
tumor growth [39,40]. Unlike the previous model, the 
immunogenic tumor model explicitly tracks the 
dynamics of the immune cells (variable x) in response to 
the tumor cells (variable y). The SDEs are as follows: 

dx{t) = (ai — a 2 x (t) + a 3 x (t)y(t)) dt 

+ (M*W -xi) + b 12 (yii) - yi)) dWi(t) 

dy{t) = (b iy (t) (1 - b 2 y(t)) - x{t)y{t)) dt 

+ (M*(0 - *i) + b 22 (y(t) - yi)) dW 2 {i) 

The parameters %\ an y 1 denote the stochastic equili- 
brium point of the model. Briefly, the model assumes 
that the amount of noise increases with the distance to 
the equilibrium point. 



For this model, we considered the following property: 
starting from 0.1 units each of tumor and immune cells, 
is there at least a 1% chance that the number of tumor 
cells could increase to 3.3 units. The property can be 
encoded into the following BLTL specification: 

Pr> 0 .oi(F 10 (x> 3.3)) 

Default model parameters were those used in [39]. Fig- 
ure 4 contrasts the number of samples needed to decide 
whether the model satisfies the property using Ltd. and 
non-LLd.. The same trends are observed as in the pre- 
vious model. That is, our algorithm requires fewer sam- 
ples than Lid. hypothesis testing, and that the difference 
between these methods increases with the Bayes factor. 

We also considered the property that the number of 
tumor cells increases to 4.0 units. We evaluated whether 
this property is true with probability at least 0.000005 
under a Bayes Factor of 100, 000. The i.i.d. sampling 
algorithm did not produce an answer even after obser- 
ving 10, 000 samples. The non-LLd. model validation 
algorithm answered affirmatively after observing 6, 736 
samples. Once again, the real impact of the proposed 
algorithm lies in uncovering rare behaviors and bound- 
ing their probability of occurrence. 



Require: stochastic differential equation Model M. AFM Specification 6. Bound on implied Radoii-Nikodym derivatives //. 
Ensure: Probability of Type I/II error is less than 6. 

Initialize n=l\ /*n is number of samples observed*/ 
Initialize :r=0; /*.r is number of successful samples*/ 
B <- 1 

while ( B > J and B < S ) do 

Choose a random implied change of probability measure Cj using an //-eventually fair testing strategy. 
Draw sample cr t} from the stochastic differential equation under the new probability measure. 

if rr n j= o then 

X n <- 1 ; x <- x + 1 
else 

X n <- 0 ; 
end if 

rl " 

j J[PiH(Xi\u)g(u) du 



Ptn(Xi\u)g(u) du 

if B > 1 then 

B ±- -\-B 

else 

B <- ij hl B 
end if 

n «— n + 1 ; 
end while 

Figure 4 Comparison of i.i.d. and non-i.i.d. sampling. Non-i.i.d. vs i.i.d. Sampling based verification for the nonlinear Immunogenic tumor 
model. 
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Discussion 

Our results confirm that non-lid. sampling reduces the 
number of samples required in the context of hypothesis 
testing - when the property under consideration is rare. 
Moreover, the benefits of non-lid. sampling increase 
with the rarity of the property, as confirmed by the 
divergence of the lines in Figures 2 and 3. We note, 
however, that if the property isn't rare then a non-lid. 
sampling strategy will actually require a larger number 
of samples than an lid. strategy. Thus, our algorithm is 
only appropriate for investigating rare behaviors. 

Conclusions 

We have introduced the first algorithm for verifying 
properties of stochastic differential equations using 
sequential hypothesis testing. Our technique combines 
Bayesian statistical model checking and non-lid. sam- 
pling and provides guarantees in terms of termination, 
and the number of samples needed to achieve those 
bounds. The method is most suitable when the behavior 
of interest is the exception and not the norm. 

The present paper only considers SDEs with indepen- 
dent Brownian noise. We believe that these results can 
be extended to handle SDEs with certain kinds of corre- 
lated noise. Another interesting direction for future 
work is the extension of these method to stochastic par- 
tial differential equations, which are used to model spa- 
tially inhomogeneous processes. Such analysis methods 
could be used, for example, to investigate properties 
concerning spatial properties of tumors, the propagation 
of electrical waves in cardiac tissue, or more generally, 
to the diffusion processes observed in nature. 
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